home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / rkqc.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-11-18  |  1.5 KB  |  68 lines
  1. /*
  2. ### quality controlled Runge-Kutta integrator one-stepper ###
  3.  
  4. Note:    a slightly modified version of rkqc in "Numerical Recipes"
  5.     c convention of [0,n-1] is used
  6. */
  7.  
  8. #include <stdio.h>
  9. #include <math.h>
  10. #define PGROW -0.20
  11. #define PSHRINK -0.25
  12. #define FCOR 6.666666666666667e-2
  13. #define SAFETY 0.9
  14. #define ERRCON 6.0e-4
  15.  
  16. void rkqc(y,dydx,n,x,htry,eps,yscal,hdid,hnext)
  17. double y[],dydx[],*x,htry,eps,yscal[],*hdid,*hnext;
  18. int n;
  19. {
  20.     int i;
  21.     double xsav,hh,h,temp,errmax;
  22.     double *dysav,*ysav,*ytemp,*dvector();
  23.     void free_dvector();
  24.     extern int model;
  25.     extern double *param;
  26.     extern int (*f_p)();
  27.  
  28.     dysav = dvector(0,n-1);
  29.     ysav = dvector(0,n-1);
  30.     ytemp = dvector(0,n-1);
  31.     xsav = (*x);
  32.     for(i=0;i<n;i++){
  33.         ysav[i] = y[i];
  34.         dysav[i] = dydx[i];
  35.     }
  36.     h = htry;
  37.     for(;;){
  38.         hh = 0.5 * h;
  39.         (void) rk4(ysav,dysav,n,xsav,hh,ytemp);
  40.         *x = xsav + hh;
  41.         (int) f_p(dydx,0,ytemp,param,*x,n);
  42.         (void) rk4(ytemp,dydx,n,*x,hh,y);
  43.         *x = xsav + h;
  44.         if(*x == xsav) system_mess_proc(1,"rkqc: Step size too small. Stop");
  45.         (void) rk4(ysav,dysav,n,xsav,h,ytemp);
  46.         errmax = 0.0;
  47.         for(i=0;i<n;i++){
  48.             ytemp[i] = y[i] - ytemp[i];
  49.             temp = fabs(ytemp[i] / yscal[i]);
  50.             if(errmax < temp) errmax = temp;
  51.         }
  52.         errmax /= eps;
  53.         if(errmax <= 1.0){
  54.             *hdid = h;
  55.             *hnext = (errmax > ERRCON ? SAFETY * h * exp(PGROW *
  56.                 log(errmax)) : 4.0 * h);
  57.             break;
  58.         }
  59.         h = SAFETY * h * exp(PSHRINK * log(errmax));
  60.     }
  61.     for(i=0;i<n;i++) y[i] += ytemp[i] * FCOR;
  62.     free_dvector(ytemp,0,n-1);
  63.     free_dvector(dysav,0,n-1);
  64.     free_dvector(ysav,0,n-1);
  65. }
  66.  
  67.  
  68.